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ABSTRACT 

We modify the public PM code developed by Anatoly Klypin and Jon Holtz- 
man to simulate cosmology with arbitrary initial power spectrum and equation 
of state of dark energy. With this tool in hand, we perform the following studies 
on the matter power spectrum. 

With an artificial sharp peak at k ~ 0.2/iMpc -1 in the initial power spectrum, 
we find that the position of the peak is not shifted by nonlinear evolution. An 
upper limit of the shift at the level of 0.02% is achieved by fitting the power 
spectrum local to the peak using a power law plus a Gaussian. This implies that, 
for any practical purpose, the baryon acoustic oscillation peaks in the matter 
power spectrum are not shifted by nonlinear evolution which would otherwise 
bias the cosmological distance estimation. We also find that the existence of a 
peak in the linear power spectrum would boost the nonlinear power at all scales 
evenly. This is contrary to what HKLM scaling relation predicts, but roughly 
consistent with that of halo model. 

We construct two dark energy models with the same linear power spectra 
today but different linear growth histories. We demonstrate that their nonlinear 
power spectra differ at the level of the maximum deviation of the corresponding 
linear power spectra in the past. Similarly, two constructed dark energy models 
with the same growth histories result in consistent nonlinear power spectra. This 
is hinting, not a proof, that linear power spectrum together with linear growth 
history uniquely determine the nonlinear power spectrum. Based on these results, 
we propose that linear growth history be included in the next generation fitting 
formulas of nonlinear power spectrum. 

For simple dark energy models parametrized by wq and w a , the existing non- 
linear power spectrum fitting formulas, which are calibrated for ACDM model, 
work reasonably well. The corrections needed are at percent level for the power 
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spe ctrum and 10% level for the derivative of the power spectrum. We find that, 
for iPeacock fe Doddsl (119961 ) fitting formula, the corrections to the derivative of 
the power spectrum are independent of w a but changing with redshift. As a short 
term solution, a fitting form could be developed for w , w a models based on this 
fact. 



Subject headings: cosmology - large-scale structure of the universe 



1. Introduction 



For cosmologists, the clustering property of matter is a gold mine. Although we are at 
the begi nning; stage of the op eration, we have already been rewarded dearly: combined with 
WMAP ( ISpergel et al.ll2003l ). SDSS data constrains cosmological parameters to reasonably 
high precision (for example, the m atter density of the universe fl m and the amplitude of the 



matt er perturbation <7s to ~ 10% (Tegma rk et al.ll2004l; lAbazaiian et al.ll2005l : ISeljak et al. 



mp 
3 



20051 )). iFosalba. Gaztanaga. k, Castanderl (120031 ) and IScranton et al.l (120031 ) found direct 



evidence of the mysterious dark energy through ISW effect using data from SDSS. Also using 
SDSS data, the baryon acoustic oscillation (BAO) feature is used to probe the expansio n 
history of the universe and cosmological distances (jEisenstein et al.l 120051 ; ICole et al.ll2005l ). 



2005 



\ made the first attempt to 



Weak lensing surveys (jJarvis. Jain. Bernstein, fc Domeyj 
constrain the time evolution of dark energy ... Moreover numerous ambitious future surveys 
of large scale structure are proposed, such as the Dark Energy S urvey (DES§),PanSTARR£0, 
Supernova/Acceleration Probe (SNAF@; Aldering et al.l 120041 ) and Large Synoptic Survey 
Telescope (LSST@). From these surveys, we want to answer questions like: Is our concordance 
cosmology model correct? What is the nature of the mysterious dark matter and dark energy? 
Do we need to modify gravity theory? Is the universe flat or curved? 



1 Presented as a dissertation to the Department of Astronomy and Astrophysics, The University of 
Chicago, in partial fulfillment of the requirements for the Ph.D. degree. 

2 Email: mzm@oddjob.uchicago.edu 

3 htt p : / / cosmology, astro . uiuc . edu/DES 

4 http: / /pan-starrs. ifa.hawaii.edu 

D http:/ /snap. lbl.gov 

6 http://www. lsst.org 
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On the other hand, we have tremendous success in theoretical modeling of structure 
formation in our universe. Starting with a handful of cosmological parame ters, we can 



comp ute the linear theory matter power spectrum with high accuracy (e.g., ISeljak et al. 



20031 ). However, we measure the nonlinear matter or galaxy power spectrum from large- 



scale structure surveys. Our understanding of the connection between linear and nonlinear 
power spectrum is improving rapidly along with our ability to simulate the complex processes 
involved. 



1.1. Baryon acoustic oscillations 

The baryon acoustic oscillation (BAO) features in the late time matter clustering, 
characterized by a single peak in the correlation function and oscillations in the matter 
power spectrum, are imprints of the acoustic oscillations in the early universe which pro- 
duces the peaks and troug hs in the CMB angular power spectrum ([Peebles k Yul Il970l ; 
Sunyaev k Zerdovichlll970l ). Although the BAO features in the matter power spectrum are 
much weaker compared to those of the CMB because the dominating dark matter did not 
participate in the acoustic oscillations, these features can be measured with high accuracy 
with sufficiently large surveys. Indeed, the fi rst detection was made in the SDSS survey in 
2005 feisenstein et alil2005l : Icole et al.lboosh . 



With the large galaxy surveys proposed, e.g. DES, PAN STARRS, SNAP, and LSST, 



and cosmological distances 


(Eisenstein, Hu, k Tegmar 


k 199S 


; Coorav. Hu, Huterer k Joffrt 




2001: Blake k Glazebrook 


2003; 


Hu k Haimanl 2003: 


Seo k Eisenstein 20031: Linder 


2003 


Matsubara 2004 : Amendola et al. 


2005: Blake k Bridle 


2005; 


Glazebrook k Blake 2005 


Dolnev et al. 


2006f). The main theoretical uncertainties, including galaxy bias, nonlinear evolution and red- 





shift s p ace distortions, are investigated in detail by IWhitd (120051). Eisenstein. Seo. k White 



(I2006f ). Eisenstein et all (120061 ) iHuff et all (120061 ) and ISmith et al.l (120061 ) . Nonlinear evolu- 



tion smears BAO peaks due to mode-mode coupling. This effect makes BAO signal weaker 
and harder to detect. An even more worrying effect is that nonlinear evolution could poten- 
tially shift the peaks which would bias the distance measurement from which the cosmological 
parameters are inferred. 



Eisenstein. Seo. k White! (120061 ) made an order of magnitude estimation of the amount 
of shift of BAO peaks due to nonlinear evolution. They use spherical collapse model to 
estimate how much the overdensity within 150 Mpc has grown if the center were over-dense; 
and the change in radius is a third of that in the over density. Their conclusion is that the 
change in scale is on the order of O(10~ 4 ). In this work, we use N-body simulations to test 
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how much the peaks are shifted due to nonlinear evolution. 



1.2. Fitting formulas of the nonlinear matter power spectrum 



Calibr ated by simulations, fitt i ng formulas of the nonlinear power spectr um have been 
developed JPeacock fc Doddsl Il996l : ISmith et aDl2003l : [McDonald et al.ll2005h . Their preci- 
sions, on the order of 10 — 20%, are sufficient for most of the current applications. Future 
surveys would have statisti cal uncertainties so low that fitting formulas with improved ac- 
curacy are highly desired (IHuterer fc Takadal 120051 ). As a step toward fulfilling this goal, 
we test the assumptions that go into the afore mentioned fitting formulas and propose the 
necessary ingredients to build a more accurate one. First, let us briefly review the existing 
fitting formulas. 



I n a series of papers ([Hamilton et al.lll99ll ; iPeacock fc Doddslll994! ; iJain. Mo. fc White 



19951 ; IPeacock fc Doddsl Il996l ). the first fitting formula of the nonlinear power spectrum is 
worked out. The foundation of the fitting formula is the so called HKLM relation which has 
two parts. The first part is to connect the nonlinear scale, r^L, to the linear scale, tl, by 
volume-averaged correlation function £(r), 



r L = [l + ^NL{r N L)} 1/3 r NL . 



The reasoning behind this is pair conservation ([Peebles Ill980l ). The second part of the HKLM 
procedure is to conjecture that the nonlinear correlation function is a universal function of 
the linear one, 

iN L {r NL ) = f NL [l L (r L )]. (2) 

If one assumes that £(r) is equivalent to matter power spectrum at an effective wavenumber 
k, then we have the power spectrum version of the HKLM relation, 



k. 



\1 



A 2 NL (k NL )f 3 k 



NL 



A 2 NL (k NL ,t) = f NL [A 2 L (k L ,t)]. 



(3) 
(4) 



It is commonly referred to as Peacock & Dodds fitting formula and f^L is calibrated using 
N-body simulation assuming ACDM cosmology. 

Another approach of ca lculating the nonlinear power spectrum is to use the halo model 
(See ISheth fc Coorayl 120021 . for a review), which splits up the contribution to the power 
spectrum into one halo and two halo terms, 



P NL (k)=P lh (k)+P m (k) 



2hi 



(5) 
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The one halo term P (k) is the contribution from two particles that reside in the same halo, 



P lh (k) = I dmn(m) 



\u(k\m)\" 



(6) 



and the two halo term is the contribution from two particles that belong to two different 
halos, 



P 2h {k) 



(771 \ 
— b(m)\u(k\m)\ 
P J 



1 2 



P L (k) . 



Here m is the mass of the halo, p is the mean de nsity of the universe, and u(k\m 



transform of the profile of a halo (NFW profile iNavarro. Frenk. fc White 



(7) 



is the Fourier 



19971. for example) 



with mass m, n(m) is the halo mass function whic h has the following form (IPress fc Schechter 
1974J : I.Tenkins et allhoOll : Isheth fc Tormenlll999l ). 



n(m) 



—f\ v )-r- 

m dm 



where v = 5 c /a(m), a(m) is the rms of the linear density field smoothed by a spherical 
top-hat filter W(k, m), 

/•oo 

dlnkA 2 L (k)W 2 (k,m) . (9) 



a 2 (m) 



Note that adding a sharp spike in the matter power spectrum would only modify a(m) below 
a mass threshold which is determined by the location of the spike. The modification to the 
mass function n{m) would also only occur below the same mass threshold. This could be 
easily shown by adding a ^-function to A|(/c) in equation[9j 



Finally, based on a fusion of the halo model and an HKLM scaling, ISmith et al.l (120031 ) 
provided an empirical fitting formula with improved perform ance (rms error bet t er tha n 7% 
is claimed). In detail, they replaced the two halo term with a lPeacock fc Doddsl (119961 ) type 
fitting form with an exponential cut off at nonlinear k to deal with translinear regime better; 
and a fitting form to the one halo term. Again, the fitting forms are calibrated using ACDM 
simulation. 

To summarize, fitting formulas of the nonlinear matter power spectrum have been de- 
veloped for ACDM cosmology with precisions on the order of 10 — 20%. The foundations 
for these fitting formulas are the HKLM scaling relation and halo model. A less obvious 
assumption behind all these fitting formulas is that the nonlinear power spectrum is deter- 
mined by the linear power spectrum at the same epoch, i.e. the history of the linear power 
spectrum does not affect the nonlinear one. 



In order to achieve the 1% accuracy required by future surveys ( IHuterer fc Takadall2005l ) . 
and with more general parametrization of dark energy models, what approach shall we take 
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to build new fitting formula of the nonlinear power spectrum? To address the issue we will 
test the assumptions that go into the existing fitting formulas using carefully constructed 
numerical experiments. In the end we will decide which assumptions should be abandoned, 
kept, and further included. 

Although all the existing fitting formulas are calibrated for ACDM cosmology, it has 
been a common practice to apply them to dark energy models in general. We will quantify 
how well this procedure works. 

The rest of the paper is organized as follows: In section 2 we briefly introduce the 
PM simulations. We demonstrate the behavior of the peak in the initial power spectrum 
in section 3, and discuss what we could learn from it. We test basic building blocks of the 
fitting formulas in section 4. In section 5 we evaluate the performance of the fitting formulas 
when applied to general dark energy models. Finally, we summarize our findings. 



2. The PM simulation 

The public PM code developed by Anatoly Klypin and Jon Holtzman is used for this 
work. The source code and manual are available at 
http : / / charon . nmsu . edu/ ~ aklypin / PM/ pmcode / 
pmcode.html 

The code assumes ACDM cosmology. In order to simulate cosmology with an arbitrary 
equation of state of dark energy, w(z), or/and arbitrary initial power spectrum, we have 
modified the code appropriately. In this section, we lay out the basic formalism used in the 
simulation and the tests we have performed. 



2.1. Formalism 

We will use r, t and v as position, time and velocity variables in physical space and x, 
t and u as the corresponding variables in comoving space. 

In physical space the Poisson equation is, 

V'$ = 47rGp, (10) 

where $ is the gravitational potential, G is the gravitational constant, and p is the density 
of the universe including all constituencies. If we define a new potential <fi as, 

= $ + i-x 2 , (11) 
2 a 
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where over-dot represents derivative with respect to r and a is the scale factor, then, in 
comoving coordinates, the Poisson equation becomes, 



= AnGpa 2 + 3- . (12) 
a 

Using Fredmann equation to substitute the a/ a term, we have, 

= 4vrGa 2 (p - p) = AitGa 2 5p , (13) 

here p is the mean density of the universe. The corresponding equations of motion in 
comoving space are, 

p = mu = mx ; p = mil = —m'V x <f) , (14) 

with m the mass of the particle. 

The PM code solves the Poisson equation [13] and the equations of motion [TH It is 
convenient to make the variables dimensionless by choosing suitable units. Anatoly Klypin's 
PM code defines the code variables (denoted by ~) as, 

x = xx , p/m=-p , <f) = (f)(f)o, p= -^Tpo , (15) 
a a 6 

with the units, 

%o = L ho J Ng /3 , p = x H , 

0o = Or o #o) 2 , p = ^$;n m . (16) 

07TG 

Equation[13] and [TU in terms of these dimensionless variables, become, 

V 2 0=|^(p-1), (17) 
2 a 



dp ri/ . ~ 7 dx . p 

/ = — = F{a)^ 

da da a' 



where 



F(a) = H /a (19) 



Q, m + Q\a 3 exp [—3 J a a =1 (l + w(a'))dlna ; ] 

Here w(a) is the equation of state of dark energy which has been generalized to an arbitrary 
function of a. 
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2.2. Brief summary of the PM code 

PM code solves equation[T7] and [18] in four main steps. 

» setup the initial positions and velocities of the dark matter particles 

» solve the Poisson equation using the density field estimated with current particle po- 
sitions 

» Advance velocities using the new potential 
» Update particle positions using new velocities 



The initial condition is setup using Zel'dovich approximation (IZel'dovichlllQTOHCrocce et al 



20061 ). The density assignment is done using Cloud- in-Cell (CIC) interpolation, which is the 
method that linearly interpolates particle mass to its eight neighboring cells. The Poisson 
equation is solved using the Fast Fourier Transform (FFT) and the particle positions and 
velocities are solved using the leapfrog integration scheme. For further details of the imple- 
mentation of these schemes, please check out the excellent write-up by Andrey Kravtsov at 



http:/ /astro. uchicago.edu/ ~andrey/Talks/PM/pm.pdf 



2.3. Code testing 

We setup a few test runs to determine the parameters to use in the simulations. These 
include the number of particles N p , number of grid cells N g , the starting redshift z- m [, evolving 
time step da and the box size L hox . 

Ideally, the larger N p and N g are the better. In practice, they are limited by the 
computation resources. For a standard run used in this work, N p = 256 3 and N g = 512 3 , 
one time step including the power spectrum evaluation takes about 2.5 min to 6min on an 
easily accessible Pentium box. Approximately 500Mb of memory is required. For a bigger 
run with iVp = 512 3 and N g = 1024 3 , which is mainly used to check the uncertainties of the 
standard runs, one step takes about 8 times longer and memory usage is on the level of 4Gb. 

Figured] shows the differences of the matter power spectrum at z = due to different 
initial redshift of the simulation, z ini . Simulations with lower z ini under estimate the nonlinear 
power spectrum. Relative to z ini = 100, z ini = 50 and z ini = 30 have p ower deficits o f 1.5% 



and 3% respectively. These are in good agreement with the results from lCrocce et al.l (120061 ). 
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Fig. 1. — The effect of initial redshift of the simulation on the matter power spectrum at 
z = 0. The initial redshift z ini and the step size da are listed in the legend. The Nyquist 
frequency is k Ny = 4/iMpc~ 1 . All of our simulations use z ini = 50 unless stated otherwise. 

For most of our applications, we use the ratio of the power spectrum. Percent level power 
deficit is not a concern for us. So we choose z- in \ = 50 as our standard initial redshift unless 
stated otherwise. As to the step size, figured shows that, with z- in i = 50, the power spectrum 
at z = differs by less than 0.3% between da = 0.004 and 0.001. Choosing da = 0.004 as 
our standard step size, a full standard run takes anywhere between 10 hours and a day. 

The box size of the simulation is more determined by the k range that a particular 
application requires than anything else. The Nyquist frequency k^ y = N g 7i/L\, ox is the 
highest mode that is resolved by the simulation. The lowest mode is 27r/L-\ 30x . Because of 
aliasing, not all the k modes in this range are faithfully calculated. To determine the precision 
of the power spectrum, we use these from simulations with higher resolution (N p = 512 3 and 
iV g = 1024 3 ) with which we compare the results from our standard simulation (N p = 256 3 
and iVg = 512 3 ). With Lt, ox = 400/?r 1 Mpc, figure[3] shows that our standard run agrees with 
the higher resolution run perfectly for k < 0.6/iMpc -1 and the agreement worsens to about 
10% at k = l.O/iMpc -1 . 

For all of our applications, except for the test of the shifting of peaks in the power 
spectrum, we use the ratio of power spectra instead of the absolute power. This practice 
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Fig. 2. — Effect of the step size of the simulation on the matter power spectrum at z — 0. 
The initial redshift z [ni = 50 and the stepsizes da are listed in the legend. The Nyquist 
frequency is k^ y = 4 /iMpc -1 . All of our simulations use da = 0.004 unless stated otherwise. 

takes care of the sample variance. FigureH] shows the uncertainty levels of the ratio. Again, 
we compare our standard run with the one that has higher resolution (N p = 512 3 and 
N g = 1024 3 ). For each simulation setup, we run two cosmological models and take the 
ratio of their power spectra. One cosmological model is ACDM (wq = w(z — 0) = — 1 and 
w a = — dw/da\ a= i = 0) and the other is w = —1 and w a = 0.3 model. According to figureHJ 
the ratios of the power spectra in simulations of different resolutions agree to better than 
0.5% all the way to k ~ 3 hMpc' 1 . 

3. the Fate of a Peak in the Initial Power Spectrum 

3.1. Setting up the experiment 

To address the issue of whether nonlinear evolution of the power spectrum would shift 
the position of BAO peaks, we construct the following numerical experiment. We put an 
artificial sharp peak in the otherwise smooth initial matter power spectrum and evolve 




Fig. 3. — Comparison o f the matter power s pectr um, at z = 0, from t he sim ulation with these 



from fitting formulas of ISmith et al.l (120031 ) and lPeacock fc Doddd ( 119961 ) . The upper panel 
shows the matter power spectra P(k) and the lower panel shows the fractional difference of 
P(k) between the higher resolution simulation with N p = 512 3 and the rest. The simulation 
box is 400 h~ 1 Mpc and the corresponding Nyquist frequencies are k^ y = 4 and 8 /iMpc -1 for 
simulations with N p = 256 3 and N p = 512 3 respectively. The initial power spectrum without 
BAO peaks is used for the Peacock & Dodds fitting formula. 



it using a PM simulation (IShandarin fc Melottl (jl990l ); iBagla fc Padmanabhanl (119971 ) used 
similar setups to test the transfer of power among different scales). The initial peak is located 
at k ~ 0.2/iMpc -1 with a FWHM = 0.013/iMpc" 1 and the maximum at k = 0.2051/iMpc -1 . 
At the maximum, P(k) is ten times what it would be without the peak. The location of 
the peak is at about where the third peak of the BAO sits and deep in the nonlinear regime 
today. The reason to use a sharp peak is to make it easier to detect peak shifting without 
using any sophisticated analysis. To set the peak amplitude high is to overcome the sample 
variance of the simulation which is at least at the level of the BAO peaks (10%). Techniques 
have been developed to get around sample variance and mak e BAO features easil y seen in 
simulations. For example, the Millennium Simulation group (ISpringel et al.l 120051 ) corrects 
the power spectrum of the actual realization of the initial fluctuations in their simulation to 
the expected input power and apply these scaling factors at all other times. It is not clear 
that these scaling factors are not affected by nonlinear evolution. To test the physics behind 
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Fig. 4. — The ratio of the matter power spectra of w a = 0.3 model to that of ACDM model 
at z = 0. Different curves are from simulations with different resolutions which are labeled 
in the legend. The Nyquist frequencies are k^ y = 4 and 8 /iMpc -1 for simulations with 
N p = 256 3 and N p = 512 3 respectively. 

it, our approach is much cleaner. 

The requirement for the simulation is such that k = 0.2/iMpc" 1 is well resolved and the 
peak is sampled with at least a few points in k space. With these considerations, we use a 
simulation with box size of 2/i _1 Gpc, number of particles iVp = 512 3 and number of grid cells 
iVg = 1024 3 . The corresponding Nyquist frequency is 1.6/iMpc" 1 and the peak is in the well 
resolved scale. With this simulation setup, the peak is sampled with five points in k space. 



3.2. The shifting of the peak 

Figure^] shows the position of the peak at z = 30 and z = 0. Tabled] lists the positions 
k and the corresponding P(k) for the points within the peak. 

The most conservative estimate of how much does the peak shift is to assume that at 
z = the maximum of the peak is somewhere between the second left (k = 0.198 /iMpc -1 ) 
and second right (k = 0.211 ft-Mpc -1 ) points relative to the highest point of the peak in 
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Fig. 5. — Testing the shift of peaks in the initial power spectrum. The blue asterisks show 
the linear power spectrum at z = 0. The nonlinear power spectrum at z = is shown using 
red plus signs. While the peak is being smeared by nonlinear evolution, the position of the 
peak is not shifted at all. 

figured We also know that the maximum of the input peak is at k = 0.204 /iMpc -1 ; So if 
the peak has shifted, the shift is less than (0.211 - 0.204)/0.204 = 3.2%. 

One reasonable assumption we can make is that there is only one peak. This narrows 
down the maximum to somewhere between the first left (k = 0.201 /iMpc -1 ) and first right 
(k = 0.207 /iMpc -1 ) points. With the input maximum known, the shift of the peak is less 
than 1.6%. 

We can certainly do better than these crude estimates. One simple method we use is 
to fit the power spectrum local to the peak with a power law plus a Gaussian. The fits for 
the power spectra at z = and z = 30 are shown in figured and figureEl The centers of the 
peak are found to be at 0.205092 hMpc' 1 and 0.205133 hMpc' 1 respectively. They agree to 
0.02%. This estimation is approaching the theoretical limit and for any practical purpose 
this is equivalent to no shift. 

These simple arguments have not taken errors into account. Sample variance does not 
apply here because the initial peak and the final peak is from the same realization. The 
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Fig. 6. — Fitting the z = 30 power spectrum near the peak with a power law plus a Gaussian. 
The peak position /c pea k and width o k are labeled in the legend. 

Poisson errors are on the level of 0.4%. This is too small to change the ordering of the points 
in the peak, so the first two estimations are not affected. The last estimation would have an 
uncertainty of ~ 10~ 6 which is too small to be relevant. 



3.3. Other effects that could shift the BAO peaks 



There are effects other than nonlinear collapse could shift the BAO peaks. The ex- 
periment above is dark matter simulation, but the actual tracer of BAO peaks are galaxies 
which we know is biased tracer of dar k matter. If g alaxy bias is scale depe ndent, it could 
potentially shift the peaks. As shown in lWhitd (120051 ). and lSmith et al.l (120061 ). this is indeed 
the case and the scale dependency should be carefully calibrated. 

Another worry is that BAO peaks sit on top of the smooth underline power spectrum 
whose slope changes due to no nlinear evolution might shift the peaks. As pointed out in 
Eisenstein. Seo. fc White! (120061 ) that this need not be the case. The argument is that one 
can easily make a template of the smooth underline power spectrum against which to achieve 
an unbiased measurement. 
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Fig. 7. — Fitting the z = power spectrum near the peak with a power law plus a Gaussian. 
The peak position /c pea k and width o k are labeled in the legend. 

3.4. Effects of the peak on nonlinear power spectrum 

To see the effects of the peak on nonlinear power spectrum, we compare the nonlinear 
power spectra evolved from initial power spectra with and without artificial peak. The 
smooth part of these initial power spectra are set to have the same amplitude. Note that 
this is different from requiring the same <y% which would allocate more power to the no peak 
case on scales outside the peak which would certainly make things more complicated and less 
clean. The simulation we use for this test is with N p = 256 3 , N g = 512 3 , Lb ox = 400/i~ 1 Mpc 
and the corresponding Nyquist frequency is k^ y = 4/iMpc -1 . We start the simulation at 
z = 30 and the step size is da = 0.004. The higher Nyquist frequency is the main reason 
that we use this run instead of the one used for the peak position test. 

The result is shown in figure[8j We see that the peak is washed out as the system 
evolves and the nonlinear power is boosted at all scales smaller than the scale where the 
peak is located. The amount of boost in power shows no scale dependency. 



The m apping between linear and nonlinear scale proposed by lHamilton et al.l (119911 ) and 
adopted by lPeacock fc Doddsl (119961 ) would have predicted that the peak being mapped to a 
nonlinear scale k^i which increases as the nonlinear power grows. We do not see that kind of 



0.1 1.0 
k [h/Mpc] 



Fig. 8. — Comparison of power spectra in simulations with or without a peak in the initial 
power spectrum. The comparison is done at a few re dshifts to show the e volut ion of the peak 
and the nonlinear power spectrum. It is clear that iPeacock fc Doddsl (I1996I ) like mapping 
between linear and nonlinear k does not exist. Instead, the effect of the peak is close to what 
halo model would predict. 



mapping here. On the other hand, the result is similar to what halo model's approach would 
predict. In the halo model, the one halo term (see equation[6]) is determined by the halo 
profile and the halo mass function n(m). The halo profile has nothing to do with whether 
there is a peak in the matter power spectrum or not. An extra peak in the linear power 
spectrum would increase the halo mass function below certain mass threshold determined by 
the position of the peak. As a result, the nonlinear power are boosted on small scales which is 
close to what we observe in this test. In detail, it does not agree with halo model's prediction 
perfectly. As shown in figure[9], halo model predicts that the fractional boost in power due 
to the existence of a peak has some weak scale dependence which is slightly different from 
what we see in the simulation. Since the addition of a spike in the power spectrum only 
changes the mass function at the low mass end where u(k\m) is scale independent at k = a 
few tenths /iMpc -1 , the one halo term acts like a shot noise term. The fractional difference 
of the power spectra is scale dependent due to the addition of this short noise like term to 
the power spectrum. We do not yet understand what accounts for this difference between 
halo model and the simulation. 
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Fig. 9. — Comparison of power spectra predicted by halo model for the cases with or without 
a peak in the linear power spectrum. The comparison is done at z = 0. The effect of the peak 
on the nonlinear power spectrum is similar to that from simulation (figure[8]), but different 
in details. 

4. Linear Growth Function and non-linear Matter Power Spectrum 



Existing fitting formulas of matter power spectrum, iPeacock fc Doddsl (119961 ) , ISmith et al. 



(120031 ) or a simple Halo model, have the following assumption: the nonlinear power spec- 
trum does not depend on the time evolution (history) of the linear power spectrum; it is 
determined by the linear power spectrum at the same epoch. 

There is a weaker version of the assumption, which states that the linear power spectrum 
plus its time evolution determine the nonlinear power spectrum. A halo model incorporating 
merger history is a good example of building fitting formula of nonlinear matter power 
spectrum under this weaker assumption. 



In this section, we test both versions of the above assumption using carefully chosen 
dark energy models. 
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4.1. Dark energy models construction 

We choose the set of binned growth suppression Gi = D + (ai)/ai and equation of state 
of dark energy Wi as parameters to be considered. Gi and Wi share the same binning with 
Az = 0.1. The linear response of Gi to the perturbations of Wi can be described by, 

5G t = Rijdwj . (20) 

Operationally the perturbation in Wi is achieved by modulating the amplitude of the form, 



S w oc exp 



2a? 



z 



(21) 



where a = 8 is used to make the shapelet sharp, and a z = 0.05 is used (2a z = 0.1 is the 
binning size). The resulting Wi binnings are slightly overlapping which does not affect the 
applications considered in this work. Gi is simply chosen as the growth suppression at the 
center of t he redshift bin zy_ The response matrix Rij is calculated numerically using the 



formula in lHu fc Eisensteinl (1l999h 



To test the strong assumption, we want to construct a w model that leads to the same 
growth factor as the w = w const model at z = but different growth in the past. In order to 
keep the linear growth factor at z = the same, the w model we are looking for should have 
an equation of state greater than w const during some epochs and smaller during some other 
epochs (to offset the effect on linear growth when it is greater). So we choose w const = —0.8, 
instead of -1, to avoid having Wi < — 1. There are many ways to achieve our goal. The 
response matrix Rij is invertible. So one could either specify the linear growth and derive 
the corresponding Wi or the other way around. Here we use the latter. First we know that 
the w model being sought after satisfies the constrain 5G\ = = RijSwj, where we explicitly 
label the redshift bin at z = as bin number 1. Furthermore, we restrict the variation of 
w to below redshift one where the effect of dark energy dominates. Finally we demand that 
w(z) > —1. We then assign values to Swj satisfying the above conditions. One of our choices 
of Wi is shown in figure[TD] (model-I). FigurefTTI shows the comparison of the corresponding 
linear growth suppression to that of w = —0.8 model. Our construction is successful since 
the two dark energy models have the same linear growth at z = and a maximum deviation 
from each other by 0.8% at z ~ 0.7. 

To test the weaker version of the assumption, we want to construct a w model that 
predicts the same (to certain precision) linear growth function as that of w = —0.8 model 
at all redshift. We apply singular value decomposition (SVD) to the response matrix R^ 
and select the eigenvectors that have small eigenvalues. These eigenvectors correspond to 
the combinations of 5wi to which 5Gi has the least response. Among the selected modes we 
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Fig. 10. — The time evolution of the equation of state of dark energy. Model-I & II have the 
same linear growth functions at z = as that of w = —0.8 model but different in the past. 
Their growth function comparisons are shown in figurefTT] Model-Ill has the same linear 
growth function (precise to 0.05% as shown in figure[T3|) . at all redshift, as that of w = —0.8 
model. 

pick the ones that involve variations of Wi at low redshift where the effect of dark energy 
dominates. Only one mode passed the selection criteria and is shown in figurefTUl as model- 
Ill. The linear growth of the selected w model agrees with that of w = —0.8 model better 
than 0.05% at all redshift as shown in figure[T51 On the other hand, the expansion histories 
of these two models differ by ~ 0.6% as shown in figurefTUl 



4.2. History does matter 

Figure[T2l and [15] summarize the results. As shown in figure[T2l the nonlinear matter 
power spectra at z = are different between w = —0.8 model and the constructed dark 
energy model-I (figure fTUl) which has the same linear growth at z = but different during 
the past. So the nonlinear power spectrum does care about the time evolution of the lin- 
ear power spectrum, contrary to what existing fitting formulas assume. Our interpretation 
is the following. The difference in the linear power spectrum during the past causes some 




-Fiducial G with w= — 0.8 

G of dark energy Model- 1 
-G of dark energy Model- II 



Fig. 11. — Comparison of the growth function of w = —0.8 model with these of the con- 
structed dark energy model-I & II which are shown in figurefTUl These dark energy models 
have the same growth functions at z — 0, but different time evolutions. The fractional 
differences shown in the lower panel are relative to w — —0.8 model. 

differences in the nonlinear power spectrum which would stay there unless some rather del- 
icate adjustments of the history of the linear power spectrum is done to offset the effect. 
An arbitrary history of the linear power spectrum does not have the ability to undo the 
difference in the nonlinear power spectrum. One possibility is that different linear growth 
historie s would produce diffe rent halo merger histories on which the halo concentration de- 



pends (IWeshsler et al.ll2005l ). So the nonlinear power spectrum which mainly comes from 
the one halo term would be different. It is hard to imagine that simply matching the linear 
growth at z = would erase the difference. So in general the difference would be there as 
shown in our test case. 

The effect of linear growth history on nonlinear power spectrum is by no means limited 
to the ~ 0.8% level shown in our test case. The nonsmooth nature of w(z) also would 
not make our argument any weaker. By a little bit of more work, we come up with another 
example that has much smoother w(z) curve (figurefTUl model-II) and matching linear growth 
at z = but maximum difference of ~ 1.5% in the past (figurefTD dash line). As shown in 
figure[T2l (dash line), the resulting nonlinear power spectra have maximum difference at ~ 2% 
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Fig. 12. — Comparison of the matter power spectrum at z = of w = —0.8 model with these 
of the constructed dark energy model-I & II which are shown in figure[TDJ These dark energy 
models have the same growth functions at z = but different time evolutions. The fractional 
differences shown in the lower panel are relative to w = —0.8 model. The simulation box 
size is Lbox = 100/iMpc -1 and the Nyquist frequency is k^ y = 16/iMpc -1 . 

level. From these examples, we see that the difference of the nonlinear power spectrum due to 
different linear growth history is at the level of the maximum deviation of the corresponding 
linear power spectra in the past. 

Our test of the weaker assumption validates it as shown in figure[15j The nonlinear 
power spectra of w — —0.8 model and the constructed dark energy model-Ill (figurefTOl 
are consistent with being essentially the same given the small difference of the linear power 
spectra. This is not a proof that the history of the linear power spectrum fully determines 
the nonlinear power spectrum, but it hints that it is possibly the case. Note that although 
there is basically no hope of distinguishing the two apparently very different dark energy 
models using their growth functions, their expansion histories differ by ~ 0.6% (see figure fTOl) 
which is close to the level that BAO could potentially have some handle on. 
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Fig. 13. — Comparison of the growth function of w = —0.8 model with that of a wiggling 
w(z) model-Ill which is shown in figurefTUl The linear growth functions of these two dark 
energy models agree to better than 0.05% at all redshift. 

5. Matter Power Spectrum and the Equation of State of Dark Energy 



With the ambitious observational efforts focused on measuring the equation of state 
of dark energy to high precision (DES,SNAP,LSST,PanSTARRS), theoretical parameter 
forecasts will need to be done more carefully than before. In particular, we have to have a 
better understanding of the nonlinear power spectrum which is very important for cosmic 
shear and any other probes that utilize the information of growth of density perturbations. 
Since numerical simulations are too expensive to run for every appli cation, fitting formulas o f 



the matter powe r spec trum are widely used. The ones prov ided by iPeacock fc Doddsl (119961 ) 



and lSmith et al.l (120031 ) are calibrated using ACDM models. McDonald et al.l (120051 ) provided 



w„ 



also 


Ma et al. 


L999h. 


Linder & White 


(2005) 



spectr um with a few percent accuracy. While the method provided b y iLinder fc White 



( 



( j2005h las not yet been wide ly adopted, a common practice is to apply IPeacock fc Doddsl 



19960 

and ISmith et al.l (120031 ) fitting formulas regardless of the dark energy models. In this 



work we evaluate how well this common practice works. 
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Fig. 14. — Comparison of the expansion history H(z)/Hq of w — —0.8 model with dark 
energy models whose equations of state are shown in figurefTUl 

Two things could go wrong if we used an inaccurate fitting formula. One is that the 
fitting formula predicts the wrong amplitude of the matter power spectrum. This is the 
part we study in this work. Another effect of incorrect fitting formula is that the predicted 
model parameter degeneracy directions and/or extent of the directions are off. Although 
not covered in this study, the latter effect is very important when the parameters are highly 



degen erate. During our study of the effect of photo-z on weak lensing tomography (IMa et al. 



20061 ) , we noticed that with no tomography binning the extent of the degeneracy among dark 



energy parameters could differ by a factor of 4 depending on whether the fitting forms of 



Peacock fc Doddsl (119961 ) or halo model is used to calculate the matter power spectrum. 



For this study we further restrict ourself to two parameter dark energy models which 
are parametrized by w and w a . 
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Fig. 15. — Comparison of the matter power spectrum at z = of w = —0.8 model with that 
of dark energy model-Ill which is shown in figurefTUl The linear growth functions of these 
two dark energy models agree to better than 0.05% at all redshift. We conclude that the 
two P(k) are consistent. 



5.1. The findings 

The quantity we use to evaluate the performance of the fitting formulas applied to 
w a 7^ dark energy models is, 



Q 



P{k,w) 



P(k,w 



(22) 



which is the fractional difference of the matter power spectrum between w a model, P(k,w), 
and the ACDM model, P(k,w = —1). The fractional difference of the power spectrum in 
terms of Q is AQ/(Q + 1) and that of the derivative of the power spectrum is AQ/Q = 
Q(fitting formula) /Q (simulation) — 1. Figures fToTfTBl show t he comparison of this qua ntify 
Q from the s imula tion with those from the fitting formulas of iPeacock fc Doddsl (1l996l ) and 
Smith et al.l (120031 ) for w a = 0.3, 0.2 and 0.1 respectively. 



First of all, the fitting formulas are doing a good job describing w a models. The fractional 
difference of the power spectrum is AQ/(Q + 1) ~ 1% and that of the derivative is AQ/Q ~ 
10% for k < 1.0 /iMpc -1 . In the context of Fisher matrix analysis, the Fisher matrix elements 
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Fig. 16. — Performance of the fitting formulas when applied to non-ACDM cosmology. In 
this particular case, Wq = —1.0 and w a = 0.3. The simulation box is 400 h~ 1 Mpc and the 
corresponding Nyquist frequency is k^ y = 4/iMpc -1 . The initial power spectrum without 
BAO peaks is used for the Peacock & Dodds fitting formula. 



would be 20% off at most. Again, just to remind ourself that we are only considering the 
effect of the amplitude of the power spectrum but not that of the shape of the power spectrum 
which is closely tied to the degeneracy direction of the model parameters. 



The most surprising result we find is that the fractional differences in Q from lPeacock fc Dodds 



(119961 ) fitting formula are independent of w a , as shown in figure[T9"l To turn this fact into 
a fitting formula of the matt er power spectrum for w a models, the time evolution of this 
potential correction factor to iPeacock &: Doddd (119961 ) fitting formula should also be inde- 
pendent of w a . Figurel20l shows this correction factor at z — 0.62. We see that the peak of 
the correction shifts to smaller scale (higher k) as one goes to higher redshift and, amazingly, 
it shows no dependency on w a . Although such regular behavior seems always have a good 
reason behind it, at the moment, we do not have a physical explanation. 

It is v ery tempting to us e such a fact to come up with a correction factor to the fitting 
formula of IPeacock fc Doddd (119961 ). The hump at the nonlinear scale in figured!] and the 
ones at high redshift can be fit very well using a Gaussian in k space. However, as we have 
pointed out throughout this work, the physical arguments supporting Peacock & Dodds' 
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Fig. 17. — Performance of the fitting formulas when applied to non-ACDM cosmology. In 
this particular case, Wq = —1.0 and w a = 0.2. The simulation box is 400 h~ 1 Mpc and the 
corresponding Nyquist frequency is k^ y = 4/iMpc -1 . The initial power spectrum without 
BAO peaks is used for the Peacock & Dodds fitting formula. 

fitting formula do not seem to hold: the linear and nonlinear scale mapping does not hold; 
the assumption that the history of the linear growth does not affect the nonlinear power 
spectrum is also being shown to be incorrect. Adding more corrections to Peacock & Dodds' 
fitting formula would only turn it into another "theory of epicycles". We do not think this 
is the right approach to build a good fitting formula. As a short term solution, this may be 
worth doing. We leave the readers to make the judgment. 



6. Discussion and Conclusions 

We modify the public PM code developed by Anatoly Klypin and Jon Holtzman to 
simulate cosmologies with arbitrary initial matter power spectra and the equation of state of 
dark energy. Using N-body simulation, we test various aspects of the matter power spectrum, 
including whether nonlinear evolution would shift the BAO peaks, the assumptions that go 
into constructing fitting formulas of the nonlinear power spectrum, and the precision of the 
existing fitting formulas when applied to dark energy models parametrized by wq and w a . 
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Fig. 18. — Performance of the fitting formulas when applied to non-ACDM cosmology. In 
this particular case, Wq = —1.0 and w a = 0.1. The simulation box is 400 h~ 1 Mpc and the 
corresponding Nyquist frequency is k^ y = 4/iMpc -1 . The initial power spectrum without 
BAO peaks is used for the Peacock & Dodds fitting formula. 

Any shift of the BAO peaks in the matter power spectrum would bias the cosmological 
distance inferred hence the cosmological parameters derived. To test whether nonlinear 
evolution would shift BAO peaks, we implement an artificial sharp peak in the initial power 
spectrum and evolve it using the PM code. We find that the position of the peak is not 
shifted by nonlinear evolution. An upper limit of the shift at the level of 0.02% is achieved 
by fitting the power spectrum local to the peak using a power law plus a Gaussian. This 
implies that, for any practical purpose, the baryon acoustic oscillation peaks in the matter 
power spectrum are not shifted by nonlinear evolution. There are other effects such as galaxy 
bias that could potentially shift the peaks. These should be carefully calibrated as well. 

We also find that the existence of a peak in the linear power spectrum would boost 
the nonlinear power at all scales evenly. This is contrary to what HKLM scaling relation 
predicts, but roughly consistent with that of halo model. The scale dependence is slightly 
different from halo model prediction in detail. Further study is required to understand the 
cause of the difference. 

All existing fitting formulas of the nonlinear power spectrum assume that the linear 



-28- 




Fig. 19.— The fractional differences of Q (defined in equation l22l . 

Q( fitting formula) / Q (simulation) — 1, for different w a . The power spectrum is 
evaluated at z — 0. 



power spectrum uniquely determines the nonlinear one, regardless of the linear growth his- 
tory. To test this assumption, we construct two dark energy models with the same linear 
power spectra today but different linear growth histories and compare their nonlinear power 
spectra from the simulation. We find that the resulting nonlinear power spectra differ at the 
level of the maximum deviation of the corresponding linear power spectra in the past. Simi- 
larly, two constructed dark energy models with the same growth histories result in consistent 
nonlinear power spectra. This is consistent with but not a proof of the conventional wisdom 
that together, the linear power spectrum and the linear growth history uniquely determine 
the nonlinear power spectrum. 

Next generation large-scale structure surveys need better fitting formulas than what are 
available now. With the high accuracy required, new fitting formulas should be based on 
solid foundations. Our results suggest that we should abandon HKLM scaling relation, keep 
halo model and further include linear growth history to build the next generation fitting 
formulas of nonlinear power spectrum. Note that our study only includes dark matter. To 
build fitting formulas with percent level accuracy , one has to includ e the effect of bary ons 



f lWhitd 12004 ; IZhan fc Knox 



2004 Ijing et al.l 120061 ) and substructure flHagan et all 120051 ). 
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Fig. 20.— The fractional differences of Q (defined in equation l22l) . 

Q(fitting formula) / Q (simulation) — 1, for different w a . The power spectrum is 
evaluated at z = 0.62. 



For simple dark energy models parametrized by Wq and w a , the existing nonlinear power 
spectrum fitting formulas, which are calibrated for ACDM model, work reasonably well. The 
corrections needed are at percent level for the power spectr um and 10% level for the deriva- 
tive of the power spectrum. We find surprisingly that, for iPeacock fc Doddsl (119961 ) fitting 
formula, the corrections needed to the derivative of the power spectrum are independent of 
w a but changing with redshift. As a short term solution, a fitting form could be developed 



for w , 



w, 



models based on this fact. 
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Table 1. Peak position and height. 





2nd left a 


1st left 


the peak 


1st right 


2nd right 


k d 


0.198 


0.201 


0.204 


0.207 


0.211 


P(k;z = 30) 


5.442 


17.75 


29.66 


26.11 


11.53 


P(k;z = 0) 


2884. 


4664. 


6265. 


5616. 


3632. 


Smearing b 


0.96 


0.48 


0.38 


0.39 


0.57 


iVolume 


48870 


51554 


54090 


54314 


56910 



a Relative to the highest point of the peak. 

b Smearing is P(k; z = 0)D 2 + (z = 30)/[P(Jfe; z = 30)D 2 + (z = 0)], 
where D + (z) is the linear growth function. 

c The number of cells in k space that contribute to the power spec- 
trum. It determines the Poisson noise. 

d Units: [Jfe^/iMpc -1 ; [P(k)}=h- 3 Mpc 3 . 



